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ABSTRACT 

The nature of structure formation around the particle free streaming scale is still far 
from understood. Many attempts to simulate hot, warm, and cold dark matter cos- 
mologies with a free streaming cutoff have been performed with cosmological particle- 
based simulations, but they all suffer from spurious structure formation at scales below 
their respective free streaming scales - i.e. where the physics of halo formation is most 
affected by free streaming. We perform a series of high resolution numerical simu- 
lations of different WDM models, and develop an approximate method to subtract 
artificial structures in the measured halo mass function. The corrected measurements 
are then used to construct and calibrate an extended Press-Schechter (EPS) model 
with sharp- A; window function and adequate mass assignment. The EPS model gives 
accurate predictions for the low redshift halo mass function of CDM and WDM mod- 
els, but it significantly under-predicts the halo abundance at high redshifts. By taking 
into account the ellipticity of the initial patches and connecting the characteristic filter 
scale to the smallest ellipsoidal axis, we are able to eliminate this inconsistency and 
obtain an accurate mass function over all redshifts and all dark matter particle masses 
covered by the simulations. As an additional application we use our model to predict 
the microhalo abundance of the standard ncutralino-CDM scenario and we give the 
first quantitative prediction of the mass function over the full range of scales of CDM 
structure formation. 



1 INTRODUCTION 



vents the haloes to form and the dark matter is distributed 



The nature of dark matter is one of the major mysteries of 
modern physics and a common point of research for par- 
ticle physics, astrophysics and cosmology. In the currently 
favored cold dark matter (CDM) model ( |Peebles|1982"| l, the 
dark matter particle is supposed to be a neutralino, the light- 



in a smooth background field instead (Smith & Markovic 



20111. Just above this characteristic scale, haloes form di- 



est stable particle in supersymetry (Jungman et al. 19961. 
With a mass around 100 GeV, neutralinos decouple very 
early and have extremely low thermal velocities, far too low 
to influence structure formation on scales relevant for galaxy 
formation. As a result we get the common picture of hier- 
archical collapse, where large haloes form through mergers 
of smaller ones, a process that spans the range from galaxy 
clusters down to microhaloes with masses of about the earth 
( jHofmann et al.|[200T| [Bertschinger|2006| ) . 

One possible alternative to CDM is the warm dark mat- 
ter (WDM) model, in which the dark matter particle is a 
sterile neutrino or gravitino (Bond & Szalay 1983 IDodel 



rectly through ellipsoidal collapse rather than through hier- 
arachical growth. 

Besides the CDM and WDM models, there are various 
alternative dark matter models such as collisional dark mat- 
ter, where the dark matter particles have a self interacting 
force (Vogclsberger et al. 2012 and references therein) or 



mixed dark matter, which consists of a mixture of cold and 
warm particles ( jMaccio et al.|2012| |Anderhalden et al.|2012| 
20131. All these models are indistinguishable from CDM 



at large scales and produce modified clustering below some 
characteristic scale. 

The nonlinear structure formation of a CDM universe 
without free streaming has been studied extensively and to 
high accuracy with high resolution cosmological simulations 



son fc Widrow|1994|[Colombi et al.|1996||Bode et al.|200iJ T 

These particles are much lighter and hence decouple later on, 
maintaining their thermal speed and influencing structure 
formation up to the scales of dwarf galaxies. While at large 
scales the collapse in WDM is hierarchical and identical to 
CDM, it becomes strongly suppressed below a characteristic 
mass scale, where the free streaming of the particles pre- 



I Davis et al.|1985 Springcl 2005) as well as analytical and 



semi-analytical approaches (Press & Schechter 1974 Bond 



et al. 1991 Sheth & Tormen 19991. However, as soon as 



free streaming effects are involved, both numerical simula- 
tion and analytical models fail to predict the correct halo 
abundances at low masses. Instead of a strong suppression, 
they produce large numbers of haloes at small masses. In the 
case of numerical simulations, the failure can be attributed 
to the fact that, if we represent the density field through a 
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discrete set of particles, then around each particle there will 
be a local gravitational sink. This sink can attract other par- 
ticles and can trigger the collapse of 'artificial' haloes even in 



the absence of cosmological perturbations (Wang & White 



2007). On the other hand, analytical approaches fail because 



they have a strong dependence on the adopted smoothing 



2006 



scale where the linear power drops to zero ( Bertschinger 
Schneider et al.|2012 |. Very recently there have been 



attempts to cure these problems with new simulation tech- 



niques (Hahn et al. 20121 as well as a modification of the 



extended Press-Schechter (hereafter EPS) approach ( Benson 
ct al. 2012). Whilst both of these approaches are promising, 
they have not demonstrated that they are uniquly conver- 
gent, nor do they currently reproduce the correct structure 
formation for small mass scales. 

In this work, we study the halo mass function in the 
presence of WDM free streaming, using both high resolu- 
tion numerical simulations and the EPS approach. We look 
at the effects of the artificial clumping in our simulations 
and propose a simple approximative method to subtract ar- 
tificial haloes in the mass function. After this correction the 
measured mass function exhibits the expected turnover and 
steep decrease towards small masses. In a second step we 
apply an EPS approach with adequate filtering and mass 
assignment, which recovers the downturn of the mass func- 
tion and gives a good match to the corrected measurements 
from our simulations. 

As an additional application of our EPS recipe, we pre- 
dict, for the first time, the neutralino-CDM mass function 
over the entire halo mass range, from the largest clusters 
down to the smallest earth-mass sized microhaloes. 

The paper is structured as follows: In §2 we take a gen- 
eral look at the free streaming and its effect on the linear 
power spectrum as well as the role of the late time thermal 
velocities. §3 is devoted to the numerical simulations of dif- 
ferent WDM cosmologies and the difficulty of artificial halo 
formation. In §4 we derive a model for the mass function 
with appropriate mass assignment and compare it to our 
simulations. This includes a method to correct for the ellip- 
ticity of initial patches in a spherically averaged Gaussian 
field. Finally, we apply our method to predict the mass func- 
tion of a neutralino-CDM cosmology in §5, and we conclude 
in 86. 



2 THE FREE STREAMING SCALE 

The thermal velocities of the dark matter particles have a 
direct influence on structure formation, since they tend to 
erase primordial perturbations below a certain scale. This 
scale depends on the mass of the dark matter particle as 
well as on its formation mechanism. 

Usually the effect of the free streaming is quantized by 
the length a particle travels before the primordial pertur- 
bations start to grow substantially, which happens to be 
around matter-radiation equality. This approximate calcu- 
lation leads to the free-streaming length 



v(t)dt 
«»(*) 



cdt 



+ 



v(t)dt 



(1) 



where tsR is the epoch when the dark matter par- 
ticles become non-relativistic, which occurs as soon as 
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Figure 1. Schematic representation of the Jeans mass Mj and 
the free streaming mass Mpg as they evolve through cosmic his- 
tory. In the red colored region the density perturbations are not 
growing because they are Jeans stable. The green colored region 
is Jeans unstable but the perturbations are completely wiped out 
due to the velocity free streaming. The free streaming scale cor- 
responds to the maximum value of the Jeans mass. 



7dm < rnuMC /3A:b. Here we have introduced the scale fac- 
tor a, the mass of the dark matter particle moM, and its 
characteristic temperature Tom- In the relativistic case, the 
mean peculiar velocity of the particle is simply v(t) ~ c. In 
the non-relativistic regime its momentum simply redshifts 
with the expansion: v tx a(t)~ . This leads to 



Af s w t-h^nr) 



i , 1 i * E Q 



(2) 



where ^h^nr) is the comoving size of the horizon at 4nr- Be- 
low the free streaming length Af 8 all perturbations are wiped 
out, the dark matter particles being in a smooth background 
density field instead. 

An alternative way of understanding the effects of free 
streaming can be obtained by following the critical Jeans 
mass through cosmic history. The linear evolution of a total 
matter perturbation may be expressed as 



d 2 5 ,_. . d5 



AnGpit) 



a 2 (t)k 2 



(3) 



where <5(r, t) — [p(r, t) — p(t)]/p(t) is the matter density per- 
turbation, p(t) is the background density of the Universe, 
er v is the dark matter velocity dispersion, and H(t) = a/a is 
the expansion rate. This expression holds on scales well be- 
low the horizon (and for non relativistic species). It may be 
noted that a necessary condition for growing mode solutions 
is that the right-hand-side of this equation stays positive. 
This leads one to introduce the effective Jeans mass 

1 3/2 



. r , s 4-7T 
Mj(t) = —p 



TV 

kj 



47T , . 

yPm(i) 



4Gp(t) 



(4) 



For M < Mj perturbations will be damped. Note that 
Eq. Q depends on the dark matter density p m as well as on 
the background density p, which includes all cosmic compo- 
nents. The two densities evolve differently, since dark matter 
becomes non relativistic well before matter radiation equal- 
ity. 

With the help of Eq. Q, it is now straight forward 
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to qualitatively trace the evolution of the Jeans criterion 
through cosmic history: In the very early universe the dark 
matter component is still relativistic (<r v ~ c) and the Jeans 
mass is growing. As soon as the dark matter becomes non 
relativistic, its thermal velocity dispersion cools in the Hub- 
ble flow (<t v ~ a ) and the Jeans mass becomes approxi- 
mately constant; during this still radiation-dominated era, 
p ~ a -4 while p m ~ a -3 . This is the case until about matter- 
radiation equality, where the evolution of the background 
density changes and the Jeans mass drops significantly. 

As soon as a mass scale becomes Jeans stable, all per- 
turbations below this mass scale are damped to zero. The 
free streaming scale can therefore alternatively be defined 
as the maximum value of the Jeans mass in cosmic history. 
This is usually the case around matter-radiation equality, 
where the Jeans scale happens to be of the order of the the 
free streaming scale defined in Eq. Q. 

The evolution of the Jeans mass and the corresponding 
free streaming mass is summarised in Fig. [I] where the red 
colored area represents scales with strictly no growing solu- 
tions, while the green colored region stands for scales with 
growing solutions but erased initial perturbations. While 
the Jeans and the free streaming mass are comparable until 
matter-radiation equality, they differ substantially in scale 
today. 



2.1 Power spectrum 

A detailed calculation of the effects of the free streaming on 
the power spectrum of initial perturbations can only be ob- 
tained by solving the coupled linearised Einstein-Boltzmann 
equations for all relevant species in the Universe. The result 
is usually given in terms of a transfer function, which is a 
mapping of the primordial perturbations from the end of in- 
flation to the moment when the first perturbations become 
nonlinear. 

In the case of a WDM universe several transfer functions 
have been proposed, taking into account different production 
mechanisms for the dark matter particle. In the following we 



use the formula presented in Viel et al. ( 2005 \ 

.WDM 1 1 / 2 



TwDin(fc) 



pCDM 



[l + (ak) 



2 111 -5/m 



with (i — 1.12 as well as 



a = 0.049 



HWDM 


-l.n 


!!wdm 


0.11 


■ h - 


1.22 


L keV 




_ 0.25 _ 




.07. 





(5) 



/i _1 Mpc, (6) 



which holds for a thermally produced dark matter candidate. 
A direct translation to the mass of a sterile neutrino is given 
by the fitting function 



4.43keV 



/mwDM\ 4/3 (SI 
\ IkeV 



)'( 



-1/3 



0.1225/ 



(7) 



The WDM free streaming introduces a characteristic 
scale of suppression, and it is convenient to define it to be 
the 'half-mode' scale at which the WDM transfer function 



drops to 1/2 (Schneider et al. 20121. The half-mode mass 
scale is given by 

47T 



AL 



W 5 



(8) 



and is about at the mass scale of a dwarf galaxy, depending 
on the exact WDM particle mass. 

In the case of a CDM cosmology with a neutralino 
dark matter candidate, the free streaming cutoff scale is 
much smaller. The transfer function obtained by |Green et al.| 
(20041 has the form 



T N (k) 



where 



3 V/ca 



exp 



(fe A ) (fe B ) 



(9) 



k A = 2.4 x 10 



/ m N \ 
UOOGeV/ 



1/2 



lOOGeV, 
(r kd /30MeV) 1/2 
1 + log(T kd /30MeV)/19.2 



hMpc" 1 ] 



(10) 



1/2 / Jl , \ 1/2 



fcB = 5 - 4xl ° 7 (l00GeV/ 

where Tkd is the kinetic decoupling temperature, which de- 
pends on the specific parameters of the supersymmetric 
model. Typical values for Tkd vary between 20 MeV and 
35 MeV, but more extreme values are possible (see |Green| 
et al.|2005| for more details). 

Owing to the fact that Eq. (|9| is not algebraically solv- 
able for k, an expression for the half-mode scale cannot 
be given. However, for a 100 GeV neutralino with a de- 
coupling temperature of 30 MeV the half-mode mass is 
M = 2.9 x lO _6 ft' 1 M , which roughly corresponds to the 
mass of the earth. 



2.2 The role of late time velocities 

Numerical simulations of cosmologies with a non negligible 
free streaming length are usually done by simply assuming a 
cutoff in the initial power spectra (as discussed above) with- 
out directly including thermal velocities in the numerical 
simulations. This approach has been applied for simulations 



in WDM (Bode et al. 2001 Lovell et al. 2011 Schneider 



et al.|20l2 I as well as simulations of tiny high redshift boxes 



in CDM (jDiemand et al.|2005| jlshiyama et al.|2010| |Ander-| 



halden & Diemand 2013). Neglecting the late time thermal 



velocity contribution is a good approximation as long as the 
Jeans mass is well below the mass resolution at the initial 
redshift of the simulation. 

After matter-radiation equality at z cq ~ 3200, the over- 
densities grow significantly while the Jeans mass drops with 
a rate of Mj ~ a~ 3//2 (as discussed in [Ji] and in the cor- 
responding Fig. [I]). This means that already at z ~ 600 
the Jeans mass is more than one order of magnitude and 
at z ~ 100 more than two orders of magnitude below the 
free streaming scale. For simulations with a typical starting 
redshift of z ~ 100, it is therefore impossible to see a di- 
rect effect of initial particle thermal velocities on the WDM 
mass function - the artificial clumping scale would be or- 
ders of magnitude larger than the thermal velocity Jeans 
scale at the initial epoch of the simulation. This also holds 
for halo density profiles, where the effects of thermal veloci- 
ties (transformation of halo cusps into central cores) are only 
observable if the late time velocities are artificially boosted 
( Villaescusa-Navarro fc Dalal|2011||Maccio et al.|20i2"l|Shao 
et al.|2013| )~ 
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Sim label 


L [h~ 1 Mpc] 


iVsim 


RS 


mwDM [kcV] 


M hul [h- 1 M Q } 


m p [h- 1 M Q ] 


l so ft[h x kpc] 


CDM_L256 


256 


1024 3 


345897 


oo 





1.18 x 10 9 


5.00 


CDM_L64 


64 


1024 3 


345897 


oo 





1.83 x 10 7 


1.25 


CDM.L16 


16 


1024 3 


345897 


oo 





2.86 x 10 5 


0.31 


WDM_ml.0_L256 


256 


1024 3 


345897 


1.0 


1.3 x 10 iu 


1.18 x 10 a 


5.00 


WDM_ml.0_L64 


64 


1024 3 


345897 


1.0 


1.3 x 10 10 


1.83 x 10 7 


1.25 


WDM_ml.0_L16a 


16 


1024 3 


345897 


1.0 


1.3 x 10 10 


2.86 x 10 5 


0.31 


WDM_ml.0_L16b 


16 


1024 3 


234786 


1.0 


1.3 x 10 10 


2.86 x 10 5 


0.31 


WDM_ml.0_L16c 


16 


1024 3 


123675 


1.0 


1.3 x 10 10 


2.86 x 10 5 


0.31 


WDM_m0.5_L256 


256 


1024 3 


345897 


0.5 


1.3 x 10 11 


1.18 x 10 9 


5.00 


WDM_m0.5_L64 


64 


1024 3 


345897 


0.5 


1.3 x 10 11 


1.83 x 10 7 


1.25 


WDM_m0.25_L256 


256 


1024 3 


345897 


0.25 


1.3 x 10 12 


1.18 x 10 9 


5.00 


WDM_m0.25_L64a 


64 


1024 3 


345897 


0.25 


1.3 x 10 12 


1.83 x 10 7 


1.25 


WDM_m0.25_L64b 


64 


1024 3 


234786 


0.25 


1.3 x 10 12 


1.83 x 10 7 


1.25 


WDM_m0.25_L64c 


64 


1024 3 


123675 


0.25 


1.3 x 10 12 


1.83 x 10 7 


1.25 


WDM_m0.25_L64d 


64 


1024 3 


012564 


0.25 


1.3 x 10 12 


1.83 x 10 7 


1.25 



Table 1. Numerical simulations used in this paper. Columns from left to right: simulation name; simulation box-size (L); particle number 
(iV s i m ); Random seed number (RS); mass of WDM particle (mwDM); half-mode mass-scale (M^ m ); ma ss of simulation particle s (m p ); 
comoving softening length (i so ft). The simulations with L = 256/i _1 Mpc have already been published in |Schneider et aL | |2012| . 



Benson et al. ( 2012 1 included the effect of late time ther- 



have been performed with PKDGRAV, a treecode with high or- 



mal velocities into their mass function calculation and found 
a very prominent effect, which changes the shape of the mass 
function at scales around the cutoff. Their model is based 



on work of Barkana et al. (2001 I, who included the effect of 



velocities in an isolated simulation of spherical collapse. The 
starting redshift of the Barkana et al. simulation, however, 
is at matter-radiation equality, at a time where all relevant 
perturbations are still extremely small and deep in the linear 
regime. Their spherical collapse simulation can therefore be 
understood as a simplified method to solve the linear Boltz- 
mann equation (where velocities have been included as well, 
see for example Viel et al.|2005 1 and, not surprisingly, leads 
to a cutoff at about the same scale as the cutoff in the WDM 
transfer function. Starting the spherical collapse much later 
at a redshift just before the relevant modes become non- 
linear, would strongly reduce the influence of the thermal 
velocities and would push the effect on the mass function to 
much smaller scales, orders of magnitude below the relevant 
half-mode scale of Eq. pi. 



3 SIMULATING THE WDM UNIVERSE 

In this section we present our set of WDM simulations and 
give details about the initial conditions, the gravity code, 
and the characteristics of the individual runs. In a second 
part we then discuss the issue of artificial clumping in detail 
and develop a way to deal with it for the purpose of the halo 
mass function. 



3.1 Characteristics of the simulations 

For all our simulations we used a WAMP7 cosmology with 
the parameters fi m = 0.2726, Qa = 0.7274, fi b = 0.046, 
h 



0.704, n s = 0.963, and a 8 = 0.809 (Komatsu et al 



20111. The CDM transfer function was generated with the 
CAMB code of Lewis et al. (20001. For the initial conditions 



we used the 2LPT code (Scoccimarro 1998 Crocce et al 



20061 with an initial redshift of zic = 49 for runs with 
box size L = 256 /i _1 Mpc and zic = 99 for runs with box 
sizes L = 64 ft _1 Mpc and L — 16 /i _1 Mpc. All simulations 



der multipole expansion and adaptive timestepping (Stadel 



2001) 



A summary of all simulations including some impor- 
tant physical and numerical quantities are listed in Table [I] 
All simulations have 1024 3 particles. The simulations with 
boxsize L = 256 ft _1 Mpc have already been published in 



pchneider et al. ( 2012 1, the others have been performed dur- 
ing the last year on the SuperMUC cluster in Munich and 
the CSCS cluster in Lugano. For the instructive purposes 
of resolving halo formation well below the cutoff scale, we 
focus on cosmologies somewhat warmer than the canonical 
2 keV WDM candidate. 

The halo finding was done using a friends-of-friends 
(FoF) algorithm ( Davis et al.|1 985) provided by the TV-Body 
ShorFI with the usual linking length of b = 0.2 and with no 
unbinding of haloes. 

For the smallest boxes of L = 16 /i _1 Mpc we performed 
a finite volume correction to compensate for the missing 
large-power modes. This was done in the simplest way by 



truncating the integrals present in Eqs ( 15 I and (161 at scales 



larger than the box length, taking the ratio of the truncated 
to non-truncated mass function and multiplying this ratio 



to the simulation measurements (see for example Watson 
et al. 2012|. Whilst this simple correction neglects the ef- 



fects of the discrete Fourier mode distribution and the run 
to run sample variance of each realization, the resulting cor- 
rected mass function is similar to that obtained using the 
correction technique of Reed et al. (20071. The result is an 



increase in the abundance of haloes at large mass scales, 
and, somewhat counter intuitively, a decrease in the abun- 
dance of haloes for smaller mass scales. For the larger boxes 
L = {64, 256} fe _1 Mpc, no finite volume correction has been 
implemented, since it has an negligible influence on the halo 
abundance. 



1 www-hpcc . astro .Washington. edu/tools/f of .html 



© 2012 RAS, MNRAS 000, [lL?? 



Halo Mass Function and the Free Streaming Scale 5 




Figure 2. Rcdshift zero mass functions measured in our simula- 
tions. Black is the CDM cosmology, blue, green and red are WDM 
cosmologies with particle masses of rrtwDM £ {0.25, 0.5, 1.0} keV. 
Circles correspond to simulations with L = 256 /i —1 Mpc, trian- 
gles to simulations with L = 64 h~ 1 Mpc and squares to simula- 
tions with L = 16 /i _1 Mpc (each with 1024 3 particles). The low 
mass upturn in the WDM runs is due to spurious structures. The 
dashed lines represent the usual Sheth-Torman mass function. 



3.2 Artificial haloes and resolution 

Numerical simulations of cosmologies with truncated ini- 
tial power spectrum produce artificial clumping at scales 
beyond the cutoff. In the simulation outputs these artifacts 
are most visible in cosmic filaments, where they form equally 



spaced clumps that are strongly resolution dependent ( Wang 
fc White||20CJ7} |Schneider et al.|[2072t >. They seem however 
to be present in all environments from underdense voids to 
high density clusters and therefore cannot simply be cut out 
of the simulation analysis. 

In the halo mass function the artificial clumps appear as 
a steep power-law, which leads to a prominent upturn below 
a characteristic mass scale. In Fig. [2] we plot the mass func- 
tion of all simulated boxes together with the usual Sheth- 
Torman prediction (dashed lines). The different colors corre- 
spond to different cosmologies (black: CDM, red: ttiwdm = 
1.0 keV, green: mwDM = 0.5 keV, blue: mwDM = 0.25 keV), 
the symbols denote the box size of the simulations (cir- 
cles: L — 256/i _1 Mpc, triangles: L = 64 h~ 1 Mpc; squares: 
L = 16 /i _1 Mpc). 

The figure clearly shows the dependence of the artif- 
ical upturn of the mass function on the resolution of the 
simulation. Furthermore, we note that the mass function of 
the artificial clumps display a power-law behaviour. Inter- 
estingly, the power-law index of the artificial mass function 
does not appear to depend strongly on the mass of the WDM 
particle. Instead, it becomes more negative with decreasing 
resolution. 

The presence of the artificial clumping is a very seri- 
ous problem for numerical simulations, since increasing the 
resolution becomes tremendously expensive. Wang & White| 
(20071 noticed that the effective converging resolution only 
goes as N 1/3 with increasing particle number N. However, 



as Fig.[2]shows things are even worse than that, because the 
slope of artificial upturn becomes shallower as the simula- 
tion resolution increases, which makes it more difficult to 
distinguish between the artificial and real part of the mass 
function. Also, it is not possible to simply shrink the box size 
of the simulation to increase the mass resolution. Owing to 
the flatness of the WDM mass function, small boxes simply 
do not have enough haloes to get a statistically meaningful 
result. The combination of all these effects makes it incredi- 
bly challenging to probe mass scales significantly below the 
half-mode mass scale. It is therefore crucial to attempt to 
obtain an improved theoretical understanding of how the 
mass function behaves at these scales. 

Additionally to the resolution problems mentioned 
above, the exact scale of the artificial upturn is not always 
easy to determine. It is possible that there are still phys- 
ical haloes in the artificially driven power-law part of the 
mass function. On the other hand the power-law of artificial 
clumps could extend to higher masses than what is appar- 
ent and change the shape of the mass function at scales well 
above the visible upturn. 

3.3 Artificial haloes and environment 

In order to examine the distribution of artificial haloes, we 
measure the mass function in different environments. If the 
visual impression that the artifacts lie predominantly in fil- 
aments is correct, then it should be possible to extract the 
halo mass function down to smaller mass scales by exclu- 
sively looking in very underdense environments. 

We measure what we will refer to as an "approximate 
conditional mass function" by imposing a halo isolation cri- 
teria. In practice, we choose a specific nearest neighbor iso- 
lation criterion and only retaining 'isolated' haloes whose 
nearest neighbour is further away than the distance d, de- 
fined in units of the box size L. This gives an approximate 
measure of the underdense conditional mass function. 'Non- 
isolated' halos on the other hand - the ones with at least 
one neighbour closer than d - can be used to approximate 
an overdense conditional mass function. 

In Fig.[3]we plot the approximate conditional mass func- 
tion of various overdense and underdense environments with 
varying distance d. The mass function measurements are 
normalised in a way that an integration over all mass bins 
leads to the same effective halo number - i.e. the measured 
mass function in a certain environment dn en v/rflogM is 
multiplied with the ratio N env /Ntot, where N cnv and JV to t 
are the total number of haloes in the specific environment 
and in the whole box, respectively. As expected from vi- 
sual appearances, the underdense regions are less contami- 
nated with spurious halos than the unconditional mass func- 
tion, as indicated by the mass scale of the artificial upturn. 
Overdense regions on the other hand, have more artifacts 
and the artificial upturn therefore happens at larger masses. 
The power-law slope of associated with the artificial haloes, 
seems to be independent of environment (c.f. Fig.p] and even 
the WDM model). 

The fact that the power-law of artifacts is 'exposed' up 
to larger masses in overdense regions means that it must 
extend far beyond the visible upturn even in the uncondi- 
tional mass function. By extension, it also implies that some 
of the halos in the artificial upturn must be real. In order to 
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Figure 3. Approximate conditional mass function of underdense 
and overdense environments for the m^DM = 0.25 keV model. 
Haloes in underdense (overdense) regions have no (at least one) 
neighbour within the distance d, defined in units of the box size 
L. The artificial upturn is shifted to smaller (larger) masses in 
underdense (overdense) regions. The power-law of the artificial 
haloes seems to be independent of the environment. The halo 
environment is defined via the distance to the nearest neigbour 



get a meaningful approximation of the physical mass func- 
tion, the artificial power-law should therefore be subtracted 
from the measurements. We do this by individually fitting 
a separate power-law function to the artificial part of every 
simulation and then subtract this power-law function from 
the measurements. This approach allows us to get a mean- 
ingful mass function down to the mass scale where the arti- 
ficial upturn is a pure power-law. However, it is important 
to note that this involves an extrapolation of the spurious 
halo mass function to large mass scales. 

In Fig.[4] the corrected mass function with power-law 
subtraction is plotted in solid bold symbols while the faint 
symbols corresponds to the original non-corrected mass 
function. The left panel shows the measured mass function 
of the WDM run with rawDM = 0.25 keV (blue), the mid- 
dle panel the the WDM run with thwdm = 0.5 keV (green) 
and the right panel the WDM run with ibwdm = 1.0 keV 
(red). The power laws, which are subtracted from the orig- 
inal mass function, are plotted as grey lines and the fitting 
of these lines is done over the yellow symbols. 

In the following we will only consider the corrected mea- 
surement of the mass function and use it to calibrate our 
analytical approaches. 



4 MASS FUNCTION IN A WDM UNIVERSE 



The EPS framework (|Press fc Schechter|[l974} |Bond et al. 
1991 Lacey fc Cole|1993 1 captures many important features 
of the end states of structure formation in the CDM model. 
In particular, the halo mass function can be defined as 



dn 



d log M 



1^ rilogq 
2M jy 1 dlogM 



(12) 



where f{v) is the first crossing distribution, a 2 (M) the vari- 
ance at the mass scale M, and p is the average density of 
the universe. On assuming uncorrelated random-walks and 
a collapse barrier set by the spherical collapse model, the 
excursion set model predicts ( Press fc Schechter|1974 Bond 
et al |1991|> : 



m = \ - 



-v/2 



s 2 c (t) 

a 2 (M) 



(13) 



where S c (t) = 1.686/ D(t) is the linearly extrapolated density 
for collapse in the spherical model and D(t) is the growth 
factor normalized to be unity at z — 0. An ellipsoidal col- 
lapse barrier gives 



-qu/2 



(14) 



where p = 0.3 and A = 0.3222 ( |Sheth fc Tormen|1999| . The 
third parameter q is predicted to be one in the ellipsoidal ap- 
proach, but Sheth & Tormen realised that the cluster abun- 
dance in simulations is better matched with the empirical 
value q = 0.707. 

In this framework all of the sensitivity of the mass func- 
tion to cosmology is encoded in the variance of the density 
perturbations on a given scale R. The variance can be ex- 
pressed as 



2 (R) = 



(15) 



where Pr^n is the linear theory matter power spectrum at 
z — and W is the Fourier transform of the filter function. 
Note that in the EPS framework the predictions are un- 
changed if we consider, rather than the density field growing 
with time and points collapsing when they cross a given den- 
sity threshold, the collapse barrier evolves with time and the 
field remains static. We adopt this latter convention. Conse- 
quently, it means that the cosmological information encoded 
in the growth of the power spectrum is transfereded to 8 c (t). 
Note also that the value of 5 C itself has a weak cosmology 
dependence ( Lahav et al.|199~ Eke et al.|1996 1, however in 
what follows this shall be neglected. 

In the case of perfectly cold dark matter, o 2 (R) rises 
monotonically towards smaller R and becomes infinite as 
R — > 0. On the other hand if free streaming takes place, the 
variance becomes constant for small but finite R, since the 



cutoff in the power spectrum truncates Eq. (j 1 5 1 
In order to evaluate Eq 



|12j we also require the loga- 
rithmic derivative of the variance with respect to the mass 
scale, given bjj^] 

dloga 2 _ 2 

~ 3ff2 



d log M 



^;P^k)W [k R) dW ^ 



(2tt)3 



dlog(fc-R)' 



(16) 



For scales well above the half-mode scale, the behavior of 
this expression is almost independent of the exact choice of 
the window function. However, this changes quite dramati- 
cally as one approaches the half-mode scale, whereupon the 
choice of the window function dictates the asymptotical be- 



havior of Eq. 1 16 1, and thereby determines the shape of the 



mass function at low masses. 



Here we assume a mass dependence of the form M <x R s 
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Figure 4. Correction of WDM mass functions for the effects of spurious structure formation in the simulations. Left panel: WDM with 
m WDM = 0.25 keV, middle panel: WDM with rrryvDM = 0.5 keV, right panel: WDM with rreyyDM = 1.0 keV. The CDM measurements 
have been added to every panel for comparison. The faint symbols correspond to the original mass function, the bold symbols correspond 
to the corrected mass function. The grey lines are the power-laws which are subtracted. The fitting is done over the yellow symbols. The 
WDM mass function decreases toward low masses only when this correction is applied. 



4.1 Choosing the right window function 

In the presence of particle free streaming the mass function 
strongly depends on the choice of the window function. Some 
common filter functions are: 

• Tophat window: This filter has the form of a sphere 
with radius -Rth and sharp boundaries in real space 
and therefore has an obvious mass assignment Mth = 
4iYpR^ H /3. In Fourier space this window has the form: 



W T h(v) = - [sin y — y cos y] . 

y 



y = kR T n- 



(17) 



The downside of the spherical tophat filter is that in Fourier 
space the random walk is correlated and the first crossing 



distribution cannot be recovered analytically (Bond et al. 
1991) |Maggiore fc Riotto||2010 \. 

• Gaussian window: This filter has no sharp bound- 
aries but the characteristic form of a Gaussian (with vari- 
ance Rqa) m rea l space. In Fourier space this property is 
maintained, leading to 



W GA (y) = e~ y2/2 , 



y = fci?GA- 



(18) 



A Gaussian window gives somewhat smoother results than 
a simple tophat window, but it has the drawback of not 
having a well denned mass. The most common practice to 
assign a mass is to normalise the filter to one in real space 
and to integrate over the filter volume. Multiplying the vol- 
ume with the average density then leads to a filter mass of 
Mga = (27r) 3 ^ 2 7rpi?Q A . The normalisation is however arbi- 
trary and this introduces an ambiguity into the mass assign- 
ment ( |Maggiore fc Riotto|2010 1. 

• Sharp-fc window: This filter is defined as a tophat 
sphere in Fourier space: 



Ws K (y) = e(l-y), 



y = kRsK, 



(19) 



and has the very appealing property that the steps of the 
random walk are uncorrelated for a Gaussian field. The 
drawback is its wiggly shape in real space, leading to contri- 
butions on all scales and making it difficult to find a reason- 
able mass assignment. The same procedure as for the Gaus- 
sian filter - i.e. normalizing the filter and integrating over 
the enclosed volume - leads to a divergent integral ( Maggiore 



fc Riotto||2010 1. Apart from the Msk oc Rg K proportional- 
ity, the mass assignment is therefore basically unconstrained 
and needs to be chosen by comparing to simulations (a more 
detailed discussion on the sharp-fc mass assignment is given 



in [4.31 



We now look at the effect of the different windows on 
the behavior of the mass function in a universe including free 
streaming. The filter specific mass function can be written 
as 



dn a 1 p 
dlogM Q _ 2 ~M a 



dlog a 2 a 



d log M Q 



(20) 



where the index a refers to a particular window function, 
e.g. a 6 {TH, GA, SK}, and so for instance /th(^th) is the 
first crossing distribution computed using the spherical top- 
hat window function. 

In the small scale limit below the half-mode scale the 



variance given in Eq. ( 15 1 becomes constant and the asymp- 



totical behavior of the mass function of Eq. ( 20 1 is therefore 
given by 



lim 



dn a 1 ^. <ilog<7 2 

o d log M a R a i? Q ->o d\ogM a 



(21) 



It is now straight forward to calculate the asymptotical be- 
havior of the mass function for different filter choices. In 
the case of a tophat window and equivalently of a Gaussian 
window we obtain 
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lim ,, dn °l r oc lim R^ 1 = oo, a = TH, GA (22) 

Hc^odlogMc Rc->-Q 

which means that the mass functions of these filters diverge 
for small scales. This result is in contradiction to our under- 
standing of structure formation around and below the half- 
mode mass scale. However, In the case of a sharp-A; window 
the asymptotical behavior of the mass function is given by 

dnsK 



lim 

«sk->« rflogMsK 



lim RskPl 
Bsk^o 



lim 



R 18 - 71 = 



In obtaining the last line we have used 



Pun(l/Rs 



^sk^cdm(1/-Rsk)2wdm(1/-Rsk) 



(23) 



(24) 



The asymptotical behavior of the sharp-A; mass function has 
the correct physical characteristics - it becomes strongly 
suppressed around the scales where free streaming is domi- 
nant. Owing to this and the fact that the first crossing dis- 
tribution was derived for the sharp-A; window, we will adopt 
this window function throughout the rest of this paper. How- 
ever, we are still left with the task of how we assign mass to 
this filter function. We shall discuss this issue in 94.31 



4.2 Mass function with sharp-A window 

It is now straight forward to derive an expression for the 
halo mass function based on the sharp-A window function. 
Since the logarithmic derivative of the window is 

dW SK 



dloe 



-y5 D (l - y), 



y = kRsK, 



(25) 



where 6 is the Dirac delta function, we can evaluate the 



integral in Eq. (16 1, obtaining 
dlogosK . 1 



Pun(l/RSK) 



dlogi?SK 27T 2 oJ K (i? S K 



^SK 



(26) 



Here we have used the fact that O(0) = 1/2. We can now 
implement Eq. ( 26 1 into the relation 



dnsi 



1 P r / rf log CT SK rflog-RsK (27) 

d log Msk 2 M SK JSK{ SK, d log R SK d log M SK ' V ' 

to obtain the mass function based on the sharp-A; win- 
dow function. Note that provided Msk oc i?| K the term 
dlogMsK/dlog-RsK = 3. 



With an appropriate mass assignment, Eq. (271 gives 



a very good match to the measurements of both the CDM 
and the WDM simulations. Especially the flattening and 
the turnover in the WDM mass function can be described 
accurately. A direct comparison to the simulations is done 
in £T4l 



4.3 Mass assignment 

In the previous sections we mentioned that the sharp-A; win- 
dow function has no well defined mass assigned to its filter 
scale. This intrinsic ambiguity can be exploited by choosing 
a mass assignment that yields a good agreement with simu- 
lations. Owing to the geometrical scaling of halo mass with 
radius (at fixed virial halo density), it is however a reason- 
able assumption to maintain the M oc i?g K proportionality 
and we can therefore write 



Sharp-k (Ellipsoidal MA) 
Sharp-k (Spherical MA) 
Tophat 




Figure 5. Mass function with tophat filter (dashed) as well as 
sharp-A; filter with spherical mass assignment (dotted) and ellip- 
soidal correction (solid). For the spherical mass assignment we 
used c = 2.7 and q = 1.0, for the ellipsoidal correction we used 
c c = 2.0 and q = 0.75. A tophat window function produces a qual- 
itatively incorrect mass function below the free streaming scale. 



, r 47T r „ 

Msk = yp[c-R S K] 



Mth 



(28) 



where c = Rtk/Rsk is a free constant. Lacey & Cole ( 1993 \ 
proposed the value c = (97r/2) 1//3 ~ 2.42, which can be ob- 
tained by normalising the filter to one in real space and 
integrating over the volume. One part of this integral is 
however diverging and Lacey & Cole set it to zero with- 



out any physical motivation (see Maggiorc & Riotto 2010 



for a more detailed discussion). In this paper, we choose 
c = 2.7, also without physical justification, in order to get 
an optimal match with our simulations. This value is not 



only larger than the one from Lacey & Cole ( 1993 1 but also 
slightly larger than the value c = 2.5 used by |Benson et al.] 
(| 2012[). The difference between our choice and the one of 



Benson et al. (20121 comes from the fact that we compare 



to the corrected mass function measurements (as explained 



in [ 3.3 1 while they compared to the direct measurement of 



the mass function and ignore all spurious haloes above the 
visible upturn. 

Furthermore, in order to obtain a good match with sim- 
ulations we also set q = 1 in Eq. (14 1. This means that with 



a sharp-A filter we can use the first crossing distribution, 
which naturally arises from the EPS approach with ellip- 
soidal collapse and we do not need any further empirical 
shift of the crossing barrier. Thus, somewhat suprisingly, we 
find that neither the rescaling of the first crossing distri- 



bution of Sheth & Tormen ( 1999 1 nor the rescaling of the 



critical overdensity done by Benson et al. ( 2012 1 is necessary. 
From a theoretical perspective this means that our sharp-A; 
model is competitive with the spherical tophat model, since 
the additional free parameter from the mass assignment is 
counterbalanced by the use of a more natural first crossing 
distribution. 
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Figure 6. Mass functions at redshift 1.1 (left panel), 2.4 (middle panel) and 4.4 (right panel). Again, the tophat mass function is plotted 
as dashed line and the sharp-fc mass function as dotted line while the solid line represents the mass function with ellipticity correction, 
which enables a much better match to the simulations at high redshifts. All the labels are the same as in Fig. [5] 



4.4 Comparison with simulations 



The halo mass function predictions from Eq. ( 27 1 together 



with the mass assignment of Eq. ( 28 1 can now be com- 
pared to the simulated CDM and WDM scenarios. We do 
this for the entire mass range spanned by our simulations 
M € [lO 9 , 10 15 ] hr x M Q and for redshifts z = and z = 4.4. 
The redshift z = 4.4 corresponds to the earliest output for 
which we have sufficiently good halo number statistics. 

Fig. [5] shows the mass function at z = 0, where the dif- 
ferent dark matter models are distinguished by color (black: 
CDM, blue: mwDM = 0.25keV, green: ibwdm = 0.5fceV, 
red: mwDM = 1.0 keV). The symbols represent the simula- 
tion measurements, with the shape depending on the box 
size of the simulation (circle: L — 256 /i _1 Mpc, triangles: 
L = 64 /i _1 Mpc, squares: L = 16 /i _1 Mpc). The spheri- 
cal tophat and sharp-fc mass functions are denoted by the 
dashed and dotted lines, respectively. The model represented 
by the solid line will be discussed in §4.5| 

For the CDM cosmology, both the spherical tophat and 
the sharp-fc mass functions give very similar predictions. 
There is a small difference for very high masses, where the 
sharp-fc mass function predicts slightly more haloe^] Nu- 
merical studies of the high mass end of the mass function 
give a halo abundance that lies in between the tophat and 



the sharp-fc model plotted here (Reed et al. 2007 Bhat 



tacharya et al. 2011 Watson et al. 20121. However, and 



most importantly, the figure also shows that the predictions 
for the sharp-fc mass function are in good agreement with 
the WDM data. On the other hand, the standard spherical 
tophat approach significantly overpredicts the halo abun- 
dances at low masses and has the wrong asymptotical be- 
havior (see discussion in [4.1 1. 



In Fig.[6]we plot the evolution of the mass function with 
redshift. The left, middle and right panels denote the re- 
sults at z = 1.1, 2.4 and 4.4, respectively. The line and 
symbol types are as in Fig.[5] Unfortunately, as one con- 
siders higher redshifts the sharp-fc model predictions begin 
to systematically underestimate the measured abundances. 

3 For CDM the difference at large masses disappears if a mass 
assignment of c = 2.42 a la Laccy & Cole (1993) is used. This is, 
however, not an option for WDM because it leads to an excess in 
the halo abundance around the half mode scale. 



This happens as soon as the WDM half-mode mass scale 
enters the exponential tail of the mass function. This ef- 
fect is most noticeable in the extreme WDM model (blue: 
mwDM = 0.25 keV) at z — 4.4 (right hand panel), where the 
discrepancy between model and data is nearly one order of 
magnitude. For a WDM model with more realistic particle 
mass, 772 wdm > IkeV, this effect is small and only becomes 
significant at very high redshifts. For example our most con- 
servative WDM model (red: ibwdm = 1 keV) seems to be 
reasonably well matched at z = 1.1 and z = 2.4, and by 
z = 4.4 there is only a slight under-prediction (visible in the 
absolute plot but not in the ratio plot). 

We conclude that the sharp-fc mass function, as given 



by Eq. ( 27 1 with mass assignment from Eq. ( 28 1 , seems to 



provide a reasonable match to simulation measurements at 
lower redshifts. However, the halo abundance of realistic 
WDM models with thwdm > 1 keV is likely to be system- 
atically underestimated somewhere above z > 5, when the 
half-mode mass scale enters the exponential tail of the mass 
function. 

In the next section we propose that an ellipticity cor- 
rection of the sharp-fc filter is required. On taking this into 
account we are able to obtain predictions that better de- 
scribe the data for the entire range of halo masses, redshifts 
and WDM dark matter particle masses that we consider in 
this work. 



4.5 Ellipticity correction 

One possible explanation for the break down of the sharp-fc 
mass function predictions for WDM model at high redshifts, 
could be that the initial patches of a Gaussian random field 
are ellipsoidal, whereas the filter function is spherical and 
so at best characterises the effective radius of the patch. 
The ellipticity of the patches change both with size and 
redshift, being more spherical (ellipsoidal) for large (small) 



mass scales and at higher (lower) redshift (Bardeen et al. 
|1986[ ). In CDM, the effect of the ellipticity can be folded 
into the mass assignment. In WDM however, ellipticity be- 
comes important as soon as the half-mode mass scale (or 
cutoff scale in Fourier space) is approached. At this scale 
only spherical perturbations survive, while ellipsoidal per- 
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Figure 7. Axes ratios a\/a 3 (solid) and 02/03 (dashed) at redshift (left panel) and redshift 4.4 (right panel) of the initial linear 
perturbations of the matter density field. The colors represent the different dark matter models (black: CDM; blue: m^DM = 0.25 keV; 
green: mwDM = 0.5 keV; red: mwDM = 1.0 keV) 



turbations do not form because their shortest axis lies below 
the half-mode mass scale, where no power is left. 

Correcting for the ellipticity of initial patches only af- 
fects the mass assignment and is independent of the first 
crossing distribution. The point is to take into account 
the discrepancy between a spherical filter and ellipsoidal 
patches, which is independent of the effect of ellipsoidal col- 
lapse, leading to the first crossing distribution of |Sheth fc| 
Tormen|(|1999|. 



The distribution of patch shapes for a Gaussian random 



field was established by Bardeen et al. (1986j), who found 
that the expected set of axis ratios could be described by: 



ax 
a 3 

02 
a 3 



(R a ) 



1 + 3e m (R a ) + p m {R a ) 

1 - 3em(R a ) +Pm(Ra) 

1 + 3e m (R a ) +Pm(R a ) 
1 - 2 Pm {R a ) 



(29) 
(30) 



where the ellipsoid axis ratios are defined such that 
(03 ^ 02 ^ ai) and where e m and p m are the ellipticity and 
prolateness parameters of an average ellipsoidal density per- 
turbation, both of which depend on the filter scale R a . In 
Appendix [A] we summarize some results from that work, 
which are salient for our application. In particular, the def- 
inition of ellipticity and prolateness as well as their depen- 
dence on the filter scale R a and redshift z are presented. 

In Fig. [7] we show the evolution of the expected axis 
ratios for a sharp-fc filter as described by Eqns (291 and 
( 30 1 . The left and right panels show the results at z = 
and z — 4.4, respectively, while the solid and dashed lines 
correspond to the ratios a\/a 3 and 02/03. The colors rep- 
resent the different dark matter scenarios (black: CDM; 
blue: mwDM = 0.25 keV; green: mwDM = 0.5 keV; red: 
mwDM = 1.0 keV). The patches are clearly more spherical 
at large scales and at high redshifts. Also the power spec- 
trum of the WDM model leads to more spherical patches 
close to the cutoff scale. 



For the case of the sharp-fc filter, the connection be- 
tween the average radius Rsk and the shortest ellipsoidal 
axis 03 can be obtained by comparing the volume of the 
filter to the exact volume of the ellipsoidal patch: 



R 



SK 



aia 2 a 3 



\a 3 a 3 J 



3 
a 3 



(to) 



(31) 



where the ratios cii /az and 02 / as and therefore £ depend on 



Rsk- On combining Eqns (29 1, (301 and (311 we get 

-Rsk 



a.3(Rs 



Z(Rs 



CORsk)' 

(1 + 3 Pn 



+ Pm) 2 



(1 - 3e m +p m )(l - 2p m ) 



1/6 



(32) 



(33) 



The function 03 (.Rsk) is bijective and can be easily inverted 
to obtain i?sK(a3). 

Now that we have a relation between the average size 
of the patch measured by a spherical filter and the effective 
smallest ellipsoidal axis, we can rewrite the mass assignment: 

4% . 



-p [c e RsK(a 3 )Y 



(34) 



where 03 is now used as the reference scale that connects 
the linear power spectrum to the final halo abundance. We 
now conjecture that as soon as 0,3 is below the half-mode 
mass scale, the corresponding ellipsoid is likely to be erased 
by the free streaming, while a spherical filter measuring an 
average radius Rsk, which is still above the half-mode mass 
scale, will count the ellipsoid as being existent. 

With the replacement of the relevant filter scale from 
Rsk to a 3 , the halo mass function can be constructed after 
the following recipe: 

(i) Compute the first crossing distribution /sk(^sk) as 
before, using the average patch radius to obtain v$k = 
^c/°"1k(^sk). No adaption to the ellipticity needs to be 
made here, since shear ellipticity is already encapsulated in 
the EPS approach with evolving barrier. 
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(ii) Calculate the average ellipticity e m and prolateness 
p m with respect to Rsk as described in Appendix [A] Deter- 
mine £(-Rsk) an d invert Eq. (321 to obtain -Rsk(«3)- 

(iii) Determine M e as well as d log og K /d log £13 and 
dlog M c /dlog a 3 . 

(iv) Construct the mass function as follows: 



dn S K 



dlo{ 


5 Me 


dlog 


„2 
°SK 


dloj 


?a3 


dloj 


;M e 



2 M c 



/0) 



*ilogoJ K rfloga 3 
d log a 3 d log M e 
P L in(l/a 3 ) 



d log a 3 



27T 2 (j| K (a3) 

3 dR S K 



da :i 



(35) 
(36) 
(37) 



Note that the last term d log M/d log 03 is not exactly 3 any- 
more because we have dropped the M oc 03 proportionality 



in Eq. ( 34 1 



In order to to obtain a good match with the simulations, 
we take c c = 2.0 for the mass assignment as well as q — 0.75 



for the first crossing distribution presented in Eq. 1 14 1. This 
means on comparison with the usual spherical tophat Sheth- 
Tormen approach, our model, with ellipticity correction, has 
one additional free parameter. This comes from the unde- 
termined relation between filter scale and mass inherent to 
the sharp-fc filter (see the discussion in j |4.3[ ). 

The predictions for our sharp-fc mass function with an 
ellipticity correction are plotted as the solid lines in Fig.[5] 
and Fig. [6] At redshift zero the differences between the re- 
sults from the sharp-fc models with spherical and ellipsoidal 
mass asignment are very small. At higher redshifts however, 
the corrected ellipsoidal mass asignment model matches well 
the estimates from the simulations. It shows a small under- 
prediction relative to the highest resolution simulations, 
however the finite volume correction is comparable to the de- 
viations between model and data. Furthermore uncertainties 
remain in our approxiamte scheme for correcting the mass 
function for spurious haloes. 

To summarize we can say that as long as the half-mode 
mass scale of the WDM model is situated well above the 
exponential part of the mass function, the sharp-fc model 
with spherical mass assignment works well and no ellipticity 
correction needs to be done. If the half-mode mass scale is 
situated in the exponential part of the mass function, then 
the sharp-fc model with spherical mass assignment under- 
predicts the halo abundance significantly and an ellipticity 
correction becomes necessary. The smaller the mass of the 
WDM particle or the higher the redshift of interest, the more 
the half-mode mass scale approaches the exponential part of 
the mass function, making a ellipticity correction essential. 
For a realistic WDM model with a particle mass around 
mwDM ~ 2 keV, the ellipticity correction is only important 
for z > 5. 

The obtained mass function has a self-similar shape 
for different cutoff scales. Thus, the half-mode mass scale 

defines a relation between the 



introduced in S2.1 



M^ 

cutoff in fc-space and the characteristic suppression scale in 
real space and can be used to pin down the peak in the mass 
function at redshift zero. Independently of the WDM model, 
the peak is located at Mp Ca k = 0.55A/^) M . Other interest- 
ing scales are also directly related to M^> M - for example, 
the scale where the halo abundance has dropped by a factor 
of 10 with respect to the peak - i.e. M10 — 0.053A/^d M . In 



principle, with perfect observations, the rate of decline of the 
halo mass function below M pea k (or for example the ratio 
Mio/Mpe a k) might be used to constrain the WDM particle 
mass. 



5 PREDICTING THE MASS FUNCTION OF A 
NEUTRALINO-CDM UNIVERSE 

In the standard ACDM picture of the universe, structure 
formation spans an enormous range of mass scales from the 
most massive galaxy super clusters down to about Earth- 
mass microhaloes. The reason for this huge hierarchy is the 
mass of the favored WIMP particle, the neutralino, which is 
of the order of 100 GeV and forms an extremely cold dark 
matter fluid at freeze-out with a tiny half-mode mass scale 
around one parsec. 

In a neutralino-CDM universe the effects of the free 
streaming are completely negligible for most astrophysical 
processes, which serves also to make their detection a chal- 
lenge. For example, microhaloes are not sufficiently compact 
to be good targets for microlensing, and are much too small 
to gravitationally compress baryons to trigger star formation 
(i.e. the cosmic baryon Jeans mass is orders of magnitude 
larger). Currently, the most promising means for detecting 
microhaloes is via indirect detection through neutralino an- 
nihilation and the resulting by-products. If the microhaloes, 
or at least their central cusps, survive the tidal forces of the 
Milky Way potential, and if they do not get disrupted by 
gravitational interactions with stars, then they are potential 



gamma-ray sources of the self-annihilating neutralino (Go 
crdt et al. 2007. Koushiappas 2009). However, the chances of 
having a microhalo close enough to the solar system for easy 
detection are low, and the overall enhancement of the signal 
is moderate ( Kamionkowski & Koushiappas 2008 . Schneider 
et al.|2010) . 

In the following section, we use pre-existing simulations 
to test our model for the neutralino-CDM mass function 
around the half-mode mass scale. Finally, we give a predic- 
tion for the neutralino-CDM mass function over the entire 
range of scales relevant for structure formation. 



5.1 Comparison with microhalo simulations 

We wish to compare the predictions from our model with 
results from JV-body simulations of the neutralino-CDM 
model. There are several numerical simulation studies pre- 
sented in the literature ( Diemand et al.| 2005 Ishiyama et al. 



2010 Anderhalden & Diemand 20131), but they tend to fo- 



cus on the halo density profile and substructure abundance 
instead of the halo mass function. The only measured halo 
mass function that we are aware of was reported by |Diemand| 
eTaLlpOOB]. 



Performing simulations of structure formation in this 
model is also incredibly challenging. Firstly, the half-mode 
mass scale is many orders of magnitude smaller than for the 
case of WDM, it being of the order M C dm ~ 10" 5 /i _1 M Q . 
Hence, one requires incredibly high resolution runs to resolve 
the relevant mass scale. Furthermore, owing to the resolu- 
tion requirement, very small simulation volumes must be 
adopted. This means that the initial power spectrum real- 
ized in the simulation volume is close to P oc fc -3 , and hence 
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tion. The dashed ones show the results from the spherical 
tophat Sheth-Tormen model. Black curves correspond to the 
CDM model with no cutoff and blue corresponds to a sce- 
nario with 100 GeV neutralino cutoff. Note that the ellip- 
ticity correction has no effect here, since we are far away 
from the exponential tail of the mass function, and using 
the simple spherical mass assignment with the usual param- 
eters c = 2.7 and q = 1 gives indistinguishable results. 

The figure clearly shows that our sharp- k mass function 
model agrees reasonably well with the simulation data and 
produces a turnover at about the right scale. Whilst this 
is an encouraging result, one should not over interpret its 
significane, owing to the poor statistics and the 'by eye' 
subtraction of artificial haloes. More detailed simulations 
would be necessary to see, for example, if the slightly steeper 
cutoff of the neutralino-CDM scenario, versus WDM, has a 
visible effect on the shape of the simulated mass function. 



Figure 8. Halo mass function of a 100 GeV neutralino CDM 
scenario with T kd = 28 MeV and WMAP1 cosmology at red- 
shift 26. The black dots correspond to the |Diemand et a.l.| ( |20C55| ) 
simulation, where artificial haloes have been removed 'by eye'. 
The dashed lines represent the Sheth-Tormen tophat mass func- 
tion, while the solid lines are our sharp-fc mass function model 
with ellipticity corrections and the usual parameters c e = 2.0 and 
q = 0.75 and yields a small scale decrease consistent with the sim- 
ulation (the simple model using spherical mass assignment with 
c = 2.7 and 9 = 1 gives indistinguishable results). Blue lines are 
with cutoff and black lines without. 



nonlinear structures grow on all scales in the box nearly si- 
mult aneouslyJe 1 gJSmither^ |Elahi et al.|2009| ) . 



Diemand et al. ( 2005 I tackled the problem by sim- 
ulating a cosmological box of size L = 3 ft _1 kpc. They 
then selected an inner region to re-simulate with a zoom 
simulation, and this was done for an effective box size of 
L = 60/i _1 pc. The initial power spectrum for their simula- 
tions was based on that from |Green et ak| ( |2004| ), with a 100 
GeV cutoff given by Eq. Q and with T kd = 28 MeV. The 
adopted background cosmology corresponded to WMAP1 
with n m = 0.268, n A = 0.732, h = 0.71 and cr 8 = 0.9 
( Spergel et al.|20 03). In order to control the nonlinear evo- 
lution of the box-scale modes, the volume was evolved from 
an initial start of z = 350 and the run halted at z — 26, by 
which time a significant fraction of the mass had collapsed 
to form microhaloes. 

Haloes were selected using a FoF algorithm and were vi- 
sually inspected. If a halo was considered to be artificial, it 
was flagged and excluded from the mass function estimation. 
Note that, because the accuracy of the visual rejection cri- 
teria has not been quantified, there is a degree of ambiguity 
left when comparing our predictions with the simulations. 

Fig. [8] presents the measured mass function from the 
simulations of Diemand et al.| ( |2005[ ) at z = 26 as a func- 
tion of halo mass. The black circles with error bars denote 
the estimates from the simulations: the three data points 
at the high mass end were derived from the the lower res- 
olution L — 3 /i -1 kpc simulation, while the five points at 
lower masses are measurements from the high resolution 
zoom runs. The solid lines show the predictions from our 
sharp-fe space model mass function with ellipticity correc- 



5.2 The complete CDM halo mass function 

We now give a prediction of the redshift zero mass function 
for a neutralino-CDM scenario for the WMAP7 cosmologi- 
cal parameters. At wavenumbers well below the neutralino 
free-streaming cut-off scale we adopt the transfer function 



of Eisenstein & Hu ( 1998). This is a good choice for our pur- 
poses, since, whilst its accuacy is 2-3% for large scale k, it 
asymptotically approaches the exact analytical solution on 
small scales (above the cutoff). The neutralino cutoff scale 
depends on the mass of the neutralino and the temperature 
of kinetic decoupling (c.f. Eqs (|9| and (111). We shall take 
the mass of the neutralino to be 100 GeV. The actual de- 
coupling temperature Tdk, however is triggered by collisions 
with leptons and is model dependent. Bcrtschi nger| ( |2006[ ) 
assumed T dk = 22.6 MeV, while |Green et al.| ( |2005| found 
Tdk = 33 MeV, based on slightly different assumptions. In 
the following we adopt these two values as benchmarks, but 
note that other values are possible, depending on the specific 
parameters of supersymmtery. 

Fig.[9]shows our prediction for the mass function of dark 
matter haloes for a 100 GeV neutralino scenario as a func- 
tion of halo mass. Our predictions cover the mass range from 
the most massive haloes M ~ 10 15 Mq until below the scale 
of one Earth mass at M ~ 1O _6 M0. For a better read- 
ability of the plot, we have excised the mass range between 
1O~ 3 M and 1O 12 M , where the mass function is essen- 
tially a power law. The blue lines correspond to the model 
scenario with a 100 GeV cutoff, while the black lines are 
without cutoff. The shaded region enclosed by the dashed 
lines, denotes the standard spherical tophat mass function 
predictions spanned by the two benchmark decoupling tem- 
peratures T dk = 22.6 MeV and T dk = 33 MeV. The blue 
shaded region enclosed by the solid lines denotes our sharp- 
k filter mass function predictions with ellipticity correction 
and the usual parameters c c = 2.0 and q = 0.75. Again, the 
simple spherical mass assignment with c = 2.7 and q — 1 
gives very similar results. 

The prediction of the neutralino-CDM mass function 
will be useful for estimations of the neutralino annihilation 
rate in the local universe. 
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Figure 9. Halo mass function of a 100 GeV neutralino WMAP7 
cosmology at redshift zero. Solid line: sharp-A; filter; dashed line: 
Shcth-Tormen tophat filter. The actual half-mode mass scale of 
a 100 GeV neutralino depends on the temperature of kinetic de- 
coupling which is model dependent. Here we have plotted the two 
benchmark values = 22.6 MeV and T^k = 33 which delimit 
the blue shaded regions, where the dark blue region represents 
our new mass function model. The black lines corresponds to a 
model without cutoff. 

6 CONCLUSIONS 

In this paper we have developed a simple analytical model 
for the halo mass function in the presence of particle free 
streaming. This is an important task since numerical simu- 
lations produce artificial clumping at the scale correspond- 
ing to the power spectrum cutoff (half-mode scale) and it is 
currently impossible to numerically test structure formation 
below this scale. In a WDM scenario, however, it is crucial 
to understand the physics of structure formation below the 
half mode scale, since it is precisely on these scales where 
the alternative dark matter scenario can be either ruled out 
or confirmed. 

Our method is based on the assumption that there is a 
direct mapping between the linear power spectrum and the 
final distribution of haloes, even in the presence of a power 
spectrum cutoff. In other words, this means that we assume 
the Extended Press-Schechter (EPS) approach to work for 
all cosmologies independent of the shape of the power spec- 
trum. We note that this assumption remains unverified; it 
is for example imaginable that an EPS approach based on 
linear perturbation theory is not sufficient as soon as the 
power spectrum deviates from its quasi power-law shape, 
where EPS is much better verified against cosmological sim- 
ulations. Instead, an approach with higher order perturba- 
tion could be necessary, or the EPS method could break 
down completely. On the longer term, it will be therefore 
indispensable to develop more accurate numerical methods 
capable of modeling the physics of gravitational structure 
formation on all scales without producing artificial clump- 
ing. 

In the following, we list some key findings of this paper 
and discuss their range of applicability: 

• In order to circumvent the problem of artificial clump- 
ing in our simulations we proposed an approximate method 



to subtract numerical artifacts from the halo mass function. 
The method is based on the observation that the mass func- 
tion of artificial haloes has the form of a steep power-law 
that continues well above the mass range where artifacts 
dominate structure formation. This power-law is then simply 
subtracted from the measurements, a correction that makes 
the mass function turn over around the half-mode scale and 
decrease toward small masses, just as is expected from theo- 
retical arguments. The power-law subtraction works surpris- 
ingly well at low redshift, where the corrected measurements 
of different box sizes are well-converged. At high redshift, 
however, the method seems to be less accurate as reflected 
by a poorer convergence between different box sizes, which 
appears to be caused in part because the artificial tail of the 
mass function is no longer such a clear power-law. 

• We have constructed a simple model for the halo mass 
function based on a sharp-A; filter with constant M oc i? 3 
mass assignment and a first crossing distribution from el- 
lipsoidal collapse. It has the same number of free parame- 
ters as the standard Sheth-Tormen approach and gives ac- 
curate predictions for CDM as well as WDM scenarios at 
low redshift. At high redshift however, as soon as the WDM 
half mode scale enters the exponential tail of the first cross- 
ing distribution, this simple method breaks down leading to 
an under-prediction of the halo abundance (Fig. [7] dotted 
lines). For a canonical 2 keV WDM candidate this happens 
at z > 5. 

• The breakdown of the simple model at high redshift 
comes from the fact that ellipsoidal patches are smoothed 
with a spherical window function, something that becomes 
important in WDM around the half-mode scale, where 
patches with high ellipticity cannot survive the free stream- 
ing. To improve upon this model, we include the effect of 
the patch ellipticity by taking the shortest ellipsoidal axes 
as the reference scale instead of the average radius that re- 
sults when a spherical window function is applied. With this 
correction, our model yields an accurate mass function for 
CDM and WDM at all redshifts tested by our simulations. 

• As a further application we used our model to predict 
the behavior of the neutralino-CDM mass function, where 
the half-mode scale lies at an extremely small scale and 
where the power spectrum cutoff has a somewhat differ- 
ent form. A comparison to simulations from Dicmand ct al. 
(2005) shows reasonable agreement, but more precise sim- 
ulations are necessary to test our mass function model in 
more detail. 

The viability of the sharp-fe mass function at small 
scales well below the half mode scale still needs to be tested 
against simulations. This requires however a completely new 
numerical approach, since all common numerical schemes 
produce artifacts at the relevant scales, and this consists of 
a formidable scientific challenge. 
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APPENDIX A: GEOMETRY OF PATCHES 

Patches around peaks of a Gaussian random field have an 
ellipsoidal form in the neighborhood around the peak. This 
observation was first done in the seminal paper of |Bardeen| 
et al. (19861, and it holds independently of the form of the 



underlying power spectrum. We will now summarize some 



of the results from Bardeen et al. (19861, and detail how 



the ellipsoidal axes are connected to the underlying power 
spectrum. 

Around a local peak the density can be approximated 
by a Taylor expansion 



A, 



[25(0) - d] 



(Al) 



where the constant density d defines an ellipsoidal patch 
and connects the eigenvalues A^ of the tensor = didjS 
to the semi major axes of an ellipsoid a^. The geometry of 
an ellipsoidal patch is characterized by the ellipticity and 
prolatness parameters 



2(Ai +A 2 + A 3 ) : 



."= jr, 2 M?\ . (A2) 



2(Ai + A 2 +A 3 )' 



From Eqs (All and ( A2 1 it is now straight forward to derive 



the ratios of the ellipsoidal axes 
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Oi _ /l + 3e + p 



1-Se + p' 



ai _ J 1 + 3e + p 
m ~ V i-2p ' 



(A3) 



that are used in the main text of this work. The elliptic- 



ity corrected mass function given by the Eqns (35 I, (361 



and ([37]) can therefore be uniquely determined, provided we 
know the distribution of ellipticity and prolateness in our 
cosmology. 

|Bardeen et al.| ( |1986[ ) found analytical prescriptions of 
the probability distributions P(e,p\x) as well as P(x\u) (Eq. 
7.6 and 7.5 in their paper), where x — (Ai + A2 + \3)/(J2 is 
the sum of the eigenvalues of the tensor divided by the 
second spectral moment. The spectral moments are defined 
as 



aj{R) = 



d 3 k 

(2tt) 3 



k 2j P hin {k)W z {kR), 



(A4) 



and (To is simply the variance defined in Eq. (15 1. The two 



probability distributions can be combined to obtain 
P{e,p\v)= dxP(e,p\x)P{x\u), 



(A5) 



which connects the ellipticity and prolateness parameters 
to the peak height v and therefore the underlying power 
spectrum. 

The distribution P{x\v) is sharply peaked and has a 
maximum at 



-yv + 



3(1 - 7 2 ) + (1.1 - 0.9 7 4 )e-^ 1 -^"^ 2 ' 2 
[3(1 - 7 2 ) + 0.45 + {-yv/2) 2 ] 1 / 2 + -yv/2 



(A6) 



where 7 = <J\I{oq02) (Bardeen et al 



1986 



Eq. 6.17). We 

ignore the distribution around a; m and just set P(x\v) = 
So(x — x m ), leading to to the simplified distribution 
P(e,p\u) = P(e,p\x m ). 



Bardeen et al. (19861 found that P(e,p\x ln ) is approxi- 



mately Gaussian (for large x) and has a maximum at 
1 30 



(5x1 + 6)V2 ' 



Pm = 



(5a& + 6) 2 



(A7) 



It is now straight forward to derive average axes ratios by 
substituting Eq. ( |A7[ | into Eq. (A3 1. The result of this calcu- 
lation is given in Fig.[7](see main text for more information). 



APPENDIX B: SUBTRACTION OF 
ARTIFICIAL HALOES AT HIGHER REDSHIFTS 

In the main text we discussed the power-law subtraction of 
the artificial haloes in detail, and we plotted the utilised 
power-law fit in Fig. [4] Here we will give the same plots for 



the correction of halo abundance at higher redshift (Fig. Bl I 
As a matter of fact, the higher the redshift the less the 
upturn has the shape of a single power-law. This makes it 
more difficult to do a proper fit and leaves space for ambi- 
guity in the model building. 
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Figure Bl. Correction of WDM halo abundance at higher redshifts. From top to bottom: z = 1.1, z = 2.4 and z = 4.4. From left to 
right: WDM with m^vDM = 0.25 keV, WDM with m^m = 0.5 keV and WDM with m^oM = 1.0 keV. The CDM measurements have 
been added to every panel for comparison. The faint symbols correspond to the original mass function, the bold symbols correspond to 
the corrected mass function. The grey lines are the power-laws which are used for the subtraction. The fitting is done over the yellow 
symbols. 
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